program main
   !============================================================================
   !                  MODULE DEPENDENCY
   !============================================================================
   !
   !     |-<- module frozen_flow -<- module Gibbs_1p4
   !     |
   !     |-<- module solvers -<- modules msi2d5 and msi2d9
   !     |
   !     |-<- module user -<- module coefficients
   !     |    |
   !     |    |->- gridboundary.dat
   !     |
   !     |-<- module grid -<- gridboundary.dat
   !     |
   !     |-<- module data -<- input_file_parameters
   !     |
   !     |-<- module blomax
   !     |
   !     |-<- module post-processing
   !     |    |
   !     |    |->- output data
   !     |
   !     |->- main program
   !
   !============================================================================

   use data
   use grid
   use solvers
   use user
   use blomax
   use postp
   use Gibbs_1p4
   use frozen_flow
   use equilibrium

   implicit none

   call get_numerical_solution

contains

   !============================================================================
   !
   ! List of subroutines
   !
   ! 1) get_numerical_solution
   ! 2) get_mass_flow_rate_and_thrust
   ! 3) save_backup
   ! 4) load_backup
   ! 5) get_one_species_constant
   ! 6) get_one_species_variable
   ! 7) get_frozen_flow
   ! 8) get_local_equilibrium
   !
   !============================================================================

   subroutine get_numerical_solution
      implicit none

      integer hora, min, seg
      real*8 :: ram
      character*2 aaux
      character*8 tcpuf

      ! Read input data from file
      ! Allocate and Initialize variables
      call initialize_variables

      ! Generates the boundary nodes of the domain. Also gives ig and Sg.
      call get_boundary_nodes(folder_input, folder_output, sim_id, nmbr, optm, &
         lid, reload, nx, ctd, ig, Sg, rcg, rt) ! Output: last four entries

      ! Write parameters in the output file
      call write_parameters(lid)

      ! Read boundary grid from a file and generates the grid
      call set_grid(folder_output, sim_id, nmbr, kg, nx, ny, a1, rt, x, y)

      ! Centroids of control volumes
      call get_real_centroids_xy(1, nx, ny, x, y, xp, yp) ! Output: last two entries

      ! Metrics calculation
      call get_metrics(nx, ny, x, y, xp, yp & ! Input
      , xe, ye, xen, yen, xk, yk, xke, yke, Jp & ! Output
      , Je, Jn, alphae, gamman, betae, betan)   ! Output

      ! Radius
      call get_radius(coord, nx, ny, y, yp, radius, re, rn, rp) ! Output: last four entries

      ! Wall temperature
      if (ccTw == 1) call set_wall_temperature(nx, Tw_cte, Twall)

      roe = 0.d0
      ron = 0.d0

      ! Flow initialization (1D isentropic flow)
      call  get_initial_guess(nx, ny, ig, modvis, beta, po, T0, gamma& ! Input
      , Rg, Sg, ye, yk, radius, rn & ! Input
      , M1D, p1D, T1D, u1D, p, T, u, v, ue & ! Output
      , un, Uce, Vcn, uin, vin, pin, Tin, Mw & ! Output
      , fm1D, Fd1D, Fpv1D, de, dn, ro, roe, ron& ! Output
      , au, ap, bv) ! Output

      av = au
      at = au

      ! Baldwin-Lomax module initialization
      if (modvis == 1 .and. modtur == 1) then
         call initialize_blomax(nx, ny, x, y)
      end if

      if (reload == 1) then
         call load_backup(nx, ny, ito, norm1, tcpuo, p, pin, pl, T, ro, roe, &
            ron, u, v, ue, ve, un, vn, Uce, Vcn, pmod, Ns, Yi) ! Output: all, except nx, ny, pmod, cmod and Ns
      end if

      !
      ! rid = residual file id
      !
      if (reload == 1) then
         open(rid, file = trim(folder_output) // trim(sim_id) // "-residual.dat", &
            POSITION = 'APPEND')
      else
         open(rid, file = trim(folder_output) // trim(sim_id) // "-residual.dat")
      end if

      ! Starts cpu time
      call cpu_time(tcpu1)

      call system("date")

      select case (pmod)
         case(1) ! One-species, constant properties

            call get_one_species_constant

         case(2) ! One-species, variable properties

            call get_one_species_variable

         case(3) ! Frozen flow

            call get_frozen_flow

         case (4) ! Local equilibrium flow

            call get_local_equilibrium

      end select

      close(rid)

      ! Finishes cpu time
      call cpu_time(tcpu2)

      tcpu = tcpuo + tcpu2 - tcpu1

      ! Backup file
      if (optm == 0) then
         call save_backup(nx, ny, itmax, norm1, tcpu, p, pin, pl, T, ro, roe, &
            ron, u, v, ue, ve, un, vn, Uce, Vcn, pmod, Ns, Yi) ! Output: none
      end if

      call get_mass_flow_rate_and_thrust(nx, ny, re, roe, u, Uce, fmi, fme, Fd) ! Output: last three entries

      ram = get_ram(sim_id, folder_output)

      hora = int(tcpu/3600)
      min = int((tcpu - hora*3600)/60)
      seg = int(tcpu - hora*3600 - min*60)
      write(tcpuf, "(i2)") hora
      write(aaux, "(i2)") min
      if (min < 10) aaux = "0" // adjustl(aaux)
      tcpuf = trim(tcpuf) // ":" // aaux
      write(aaux, "(i2)") seg
      if (seg < 10) aaux = "0" // adjustl(aaux)
      tcpuf = trim(tcpuf) // ":" // aaux

      write(lid,*)
      write(lid,*)
      write(lid,*) " ==============================  MAIN RESULTS  ", &
         "=============================="
      write(lid,*)
      write(lid,"(ES23.16, ' = ', A)") fmi, &
         "fmi:   Mass flow rate in the entrance (kg/s)"
      write(lid,"(ES23.16, ' = ', A)") fme, &
         "fme:   Mass flow rate in the exit (kg/s)"
      write(lid,"(ES23.16, ' = ', A)") Fd, "Fd:    Dynamic thrust (N)"
      write(lid,"(ES23.16, ' = ', A)") tcpu, "tcpu:  CPU time (s)"
      write(lid,"(ES23.16, ' = ', A)") ram, "RAM:   Memory (MB)"
      write(lid,"(ES23.16, ' = ', A)") dt, "dt1:   Final time step (s)"
      write(lid,"(I23, ' = ', A)") it-2, "it:    Number of iterations"
      write(lid,"(15x, A, ' = ', A)") tcpuf, "tcpuf: CPU time (hh:mm:ss)"

      call post_processing(nx, ny, sem_a, sem_g, w_g, w_cam, ccTw, ig, lid, &
         sim_id, x, y, rp, re, xe, ye, xk, yk, Jp, Sg, rcg, gamma, pr, go, &
         po, Rg, tcpuo, tcpu1, tcpu2, tcpu, Fpv1D, Fd1D, fm1D, fme, cp, vlp, &
         kp, gcp, Twall, u1D, p1D, T1D, M1D, uin, vin, pin, Tin, roe, Uce, &
         Vcn, de, dn, pl, bp, xp, yp, u, v, p, T, ro, pmod, Ns, species, Yi, &
         specie, optm, folder_input, folder_output, modtur, vtp)

      call get_east_boundary_field(nx, ny, u, ube)
      call get_east_boundary_field(nx, ny, v, vbe)
      call get_east_boundary_field(nx, ny, T, Tbe)
      call get_east_boundary_field(nx, ny, p, pbe)

      call get_west_boundary_field(nx, ny, u, ubw)
      call get_west_boundary_field(nx, ny, v, vbw)
      call get_west_boundary_field(nx, ny, T, Tbw)
      call get_west_boundary_field(nx, ny, p, pbw)

      call get_north_boundary_field(nx, ny, u, ubn)
      call get_north_boundary_field(nx, ny, v, vbn)
      call get_north_boundary_field(nx, ny, T, Tbn)
      call get_north_boundary_field(nx, ny, p, pbn)

      call get_south_boundary_field(nx, ny, u, ubs)
      call get_south_boundary_field(nx, ny, v, vbs)
      call get_south_boundary_field(nx, ny, T, Tbs)
      call get_south_boundary_field(nx, ny, p, pbs)

      write(15,*) it - 2
      write(15,*) tcpu
      write(15,*) u
      write(15,*) v
      write(15,*) T
      write(15,*) p
      write(15,*) ube
      write(15,*) ubw
      write(15,*) ubn
      write(15,*) ubs
      write(15,*) vbe
      write(15,*) vbw
      write(15,*) vbn
      write(15,*) vbs
      write(15,*) Tbs
      write(15,*) Tbn
      write(15,*) Tbw
      write(15,*) Tbe
      write(15,*) pbs
      write(15,*) pbn
      write(15,*) pbw
      write(15,*) pbe

      close(15)

   end subroutine get_numerical_solution

   !============================================================================

   !> Calculates the field F over the east boundary
   subroutine get_east_boundary_field(nx, ny, F, Fbe)
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension(nx*ny), intent(in)  :: F   !< A generic field
      real(8), dimension(ny),    intent(out) :: Fbe !< The field over the east boundary

      ! Inner variables

      integer :: i, j, np, npe ! Dummy variables

      i = nx - 1

      do j = 2, ny-1

         np   = nx * (j-1) + i
         npe  = np + 1

         Fbe(j) = 0.5d0 * ( F(np) + F(npe) )

      end do

   end subroutine get_east_boundary_field

   !============================================================================

   !> Calculates the field F over the west boundary
   subroutine get_west_boundary_field(nx, ny, F, Fbw)
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension(nx*ny), intent(in)  :: F   !< A generic field
      real(8), dimension(ny),    intent(out) :: Fbw !< The field over the west boundary

      ! Inner variables

      integer :: i, j, np, npw ! Dummy variables

      i = 2

      do j = 2, ny-1

         np   = nx * (j-1) + i
         npw  = np - 1

         Fbw(j) = 0.5d0 * ( F(np) + F(npw) )

      end do

   end subroutine get_west_boundary_field

   !============================================================================

   !> Calculates the field F over the north boundary
   subroutine get_north_boundary_field(nx, ny, F, Fbn)
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension(nx*ny), intent(in)  :: F   !< A generic field
      real(8), dimension(nx),    intent(out) :: Fbn !< The field over the north boundary

      ! Inner variables

      integer :: i, j, np, npn ! Dummy variables

      j = ny-1

      do i = 2, nx-1

         np   = nx * (j-1) + i
         npn  = np + nx

         Fbn(i) = 0.5d0 * ( F(np) + F(npn) )

      end do

   end subroutine get_north_boundary_field

   !============================================================================

   !> Calculates the field F over the south boundary
   subroutine get_south_boundary_field(nx, ny, F, Fbs)
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), dimension(nx*ny), intent(in)  :: F   !< A generic field
      real(8), dimension(nx),    intent(out) :: Fbs !< The field over the south boundary

      ! Inner variables

      integer :: i, j, np, nps ! Dummy variables

      j = 2

      do i = 2, nx-1

         np   = nx * (j-1) + i
         nps  = np - nx

         Fbs(i) = 0.5d0 * ( F(np) + F(nps) )

      end do

   end subroutine get_south_boundary_field

   !============================================================================

   subroutine get_mass_flow_rate_and_thrust(nx, ny, re, roe, u, Uce, fmi, fme, &
         Fd)
      implicit none
      integer, intent(in) :: nx  ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny  ! Number of volumes in eta direction (real + fictitious)
      real(8), dimension(nx*ny), intent(in) :: re   ! Radius of the center of east face of volume P
      real(8), dimension(nx*ny), intent(in) :: roe  ! Absolute density at east face
      real(8), dimension(nx*ny), intent(in) :: u    ! Cartesian velocity of the last iteraction
      real(8), dimension(nx*ny), intent(in) :: Uce  ! Contravariant velocity U at east face
      real(8), intent(out) :: fmi ! Mass flow rate at entrance
      real(8), intent(out) :: fme ! Mass flow rate at exit
      real(8), intent(out) :: Fd  ! Thrust

      real(8), parameter :: pi = acos(-1.d0)

      integer :: i, j, np

      fmi = 0.d0

      i = 1

      do j = 2, ny-1

         np = nx*(j-1) + i

         fmi = fmi + re(np)*roe(np)*Uce(np)

      end do

      fmi = fmi*2.d0*pi

      fme = 0.d0

      Fd = 0.d0

      i = nx-1

      do j = 2, ny-1

         np = nx*(j-1) + i

         fme = fme + re(np)*roe(np)*Uce(np)

         Fd = Fd + re(np)*roe(np)*Uce(np)*(u(np) + u(np + 1))/2.d0

      end do

      fme = fme*2.d0*pi

      Fd = Fd*2.d0*pi

   end subroutine get_mass_flow_rate_and_thrust

   !============================================================================

   ! Gets the current process memory in Mb
   real(8) function get_ram(sim_id, folder_output)
      implicit none
      character(len = *), intent(in) :: sim_id ! Simulation identification
      character(200) :: folder_output

      ! Inner variables
      character(len = 20) :: pid
      character(len = 200) :: fout
      real(8) :: ram

      write(pid,*) getpid()

      fout = trim(folder_output) // trim(adjustl(sim_id)) // &
         "-memory.dat"

      call system("ps -o drs -p " // trim(adjustl(pid)) // " | tail -n 1 > " &
         // trim(adjustl(fout)))

      open(10, file = fout)

      read(10,*) ram

      close(10)

      get_ram = ram/1024.d0

   end function

   !============================================================================

   ! Subroutine 03
   subroutine save_backup(nx, ny, it, norm1, tcpu, p, pin, pl, T, &
         ro, roe, ron, u, v, ue, ve, un, vn, Uce, Vcn, pmod, Ns, Yi) ! Output: none
      implicit none
      integer, intent(in) :: pmod ! Physical model
      integer, intent(in) :: Ns   ! Number of chemical species
      integer, intent(in) :: nx   ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny   ! Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: it   ! Iteraction number when backup was made
      real(8), intent(in) :: norm1 ! Norm of the linear systems of the first iteraction
      real(8), intent(in) :: tcpu  ! Elapsed cpu time
      real(8), dimension(nx*ny), intent(in) :: p   ! Pressure at center of volume P
      real(8), dimension(ny), intent(in) :: pin ! Pressure in the entrance
      real(8), dimension(nx*ny), intent(in) :: pl  ! Pressure correction at center of volume P
      real(8), dimension(nx*ny), intent(in) :: T   ! Temperature at center of volume P
      real(8), dimension(nx*ny), intent(in) :: ro  ! Specific mass at center of volume P
      real(8), dimension(nx*ny), intent(in) :: roe ! Absolute density at east face
      real(8), dimension(nx*ny), intent(in) :: ron ! Absolute density at north face
      real(8), dimension(nx*ny), intent(in) :: u   ! Cartesian velocity at center of volume P
      real(8), dimension(nx*ny), intent(in) :: v   ! Cartesian velocity at center of volume P
      real(8), dimension(nx*ny), intent(in) :: ue  ! Cartesian velocity u at center of east face
      real(8), dimension(nx*ny), intent(in) :: ve  ! Cartesian velocity v at center of east face
      real(8), dimension(nx*ny), intent(in) :: un  ! Cartesian velocity u at center of north face
      real(8), dimension(nx*ny), intent(in) :: vn  ! Cartesian velocity v at center of north face
      real(8), dimension(nx*ny), intent(in) :: Uce ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(in) :: Vcn ! Contravariant velocity V at north face
      real(8), dimension(:,:), intent(in) :: Yi  ! Mass fractions of chemical species

      integer :: i

      open(10, file = trim(folder_output) // trim(adjustl(sim_id)) // ".bkp")

      write(10,*) it, norm1, tcpu, p, pin, pl, T, ro, roe, ron, u, v, ue, ve, &
         un, vn, Uce, Vcn

      if (pmod >= 3) then
         do i = 1, Ns
            write(10,*) Yi(:, i)
         end do
      end if

      close(10)

   end subroutine save_backup

   !============================================================================

   subroutine load_backup(nx, ny, ito, norm1, tcpu, p, pin, pl, T, ro, roe, &
         ron, u, v, ue, ve, un, vn, Uce, Vcn, pmod, Ns, Yi) ! Output: all, except nx, ny, pmod, cmod and Ns
      implicit none
      integer, intent(in) :: nx   ! Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny   ! Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: pmod ! Physical model
      integer, intent(in) :: Ns   ! Number of chemical species
      integer, intent(out) :: ito ! Iteraction number when backup was made
      real(8), intent(out) :: norm1 ! Norm of the linear systems of the first iteraction
      real(8), intent(out) :: tcpu  ! Elapsed cpu time
      real(8), dimension(nx*ny), intent(out) :: p   ! Pressure at center of volume P
      real(8), dimension(ny), intent(out) :: pin ! Pressure in the entrance
      real(8), dimension(nx*ny), intent(out) :: pl  ! Pressure correction at center of volume P
      real(8), dimension(nx*ny), intent(out) :: T   ! Temperature at center of volume P
      real(8), dimension(nx*ny), intent(out) :: ro  ! Specific mass at center of volume P
      real(8), dimension(nx*ny), intent(out) :: roe ! Absolute density at east face
      real(8), dimension(nx*ny), intent(out) :: ron ! Absolute density at north face
      real(8), dimension(nx*ny), intent(out) :: u   ! Cartesian velocity at center of volume P
      real(8), dimension(nx*ny), intent(out) :: v   ! Cartesian velocity at center of volume P
      real(8), dimension(nx*ny), intent(out) :: ue  ! Cartesian velocity u at center of east face
      real(8), dimension(nx*ny), intent(out) :: ve  ! Cartesian velocity v at center of east face
      real(8), dimension(nx*ny), intent(out) :: un  ! Cartesian velocity u at center of north face
      real(8), dimension(nx*ny), intent(out) :: vn  ! Cartesian velocity v at center of north face
      real(8), dimension(nx*ny), intent(out) :: Uce ! Contravariant velocity U at east face
      real(8), dimension(nx*ny), intent(out) :: Vcn ! Contravariant velocity V at north face
      real(8), dimension(:,:), intent(out) :: Yi  ! Mass fraction of chemical species

      ! Auxiliary variables

      logical :: lexist
      integer :: i     ! Counter

      inquire(file = trim(folder_output) // trim(adjustl(sim_id)) // &
         ".bkp", exist = lexist)

      if (lexist) then

         open(10, file = trim(folder_output) // trim(adjustl(sim_id)) // ".bkp")

         read(10,*) ito, norm1, tcpu, p, pin, pl, T, ro, roe, ron, u, v, ue, &
            ve, un, vn, Uce, Vcn

         if (pmod >= 3) then
            do i = 1, Ns
               read(10,*) Yi(:, i)
            end do
         end if

         close(10)

         ito = ito + 1

      else

         write(*,*) 'load_backup: Backup file does not exist.'

      end if

   end subroutine load_backup

   !============================================================================

   subroutine get_one_species_constant

      implicit none
      integer :: im
      integer :: cont = 0
      real*8 :: dFd

      51 format(" ==========================================================", &
         "==================== ", /, " ============ ERROR: Numero nao ", &
         "definido - Tente outros parametros ============ ", /, &
         " =================================================================", &
         "============= ", /)

      do it = ito, itmax

         ! Parada com base na tolerância para o dFd*
         if (cont > ittol) then
            exit
         else
            if (dFd < tolerance) then
               cont = cont + 1
            else
               cont = 0
            end if
         endif

         norm = 0.d0

         ! Calculating mean variables used in the eddy viscosity calculation
         if (modvis == 1 .and. modtur == 1) then
            ro_av = 0.5d0*(ro + roa)
            ron_av = 0.5d0*(ron + rona)
            u_av = 0.5d0*(u + ua)
            v_av = 0.5d0*(v + va)
         end if

         ! Updating variables for the new time step
         ua = u
         va = v
         uea = ue
         vea = ve
         una = un
         vna = vn
         pa = p
         Ta = T
         roa = ro
         rona = ron
         Ucem = Uce
         Vcnm = Vcn
         pina = pin

         ! Updating conditions at the nozzle entrance
         call get_uin_vin_pin_Tin_Mw(nx, ny, gamma, Rg, po, T0, u & ! Input
         , uin, vin, pin, Tin, Mw)        ! Output

         call get_plin_and_p_fictitious(nx, ny, pina, pin, plin, p) ! Output: last two entries

         ! Calculation of the gas properties
         call set_cp_and_gamma(nx, ny, Rg, cp_cte, cp, gcp) ! Output: last two entries

         if (modvis == 1) then

            call set_laminar_viscosity_at_nodes(nx, ny, visc_cte, vlp) ! Output: last entry

            call get_laminar_viscosity_at_faces(nx, ny, vlp, vle, vln) ! Output: last two entries

            if (modtur == 1) then

               ! Calculates the dimensional eddy viscosity at center of real volumes
               call get_eddy_viscosity_blomax(nx, ny, x, y, xp, yp, xk, yk, &
                  Jp, Jn, vln, ro_av, ron_av, u_av, v_av, vtp) ! Last one is output

               ! Extrapolation of the eddy viscosity to the fictitious volumes
               call get_eddy_viscosity_fictitious(nx, ny, vtp)

               vlp = vlp + vtp

               ! CAUTION: In order to apply correctly this function, laminar viscosity vlp must be known
               ! in real and fictitious volumes (except CV of corners SW, SE, NW and NE)
               call get_laminar_viscosity_at_faces(nx, ny, vlp, vle, vln) ! Output: last two entries

            end if

            call set_thermal_conductivity_at_nodes(nx, ny, k_cte, kp)  ! Output: last entry

            if (modtur == 1) then
               kp = kp + vtp*cp/0.9d0
            end if

            call get_thermal_conductivity_at_faces(nx, ny, kp, ke, kn) ! Output: last two entries

         end if

         ! Calculating u coefficients and source
         call get_u_coefficients(nx, ny, modvis, dt, rp, re, rn, Jp, Je, Jn &
            , ye, yk, alphae, betae, betan, gamman &
            , vle, vln, roe, ron, roa, Uce, Vcn, au)    ! Output: last entry

         call get_u_source(nx, ny, modvis, beta, dt, rp, re, rn&
            , xe, ye, xk, yk, xke, yke, xen, yen &
            , Jp, Je, Jn, roe, ron, roa, p, vle, vln &
            , Uce, Vcn, ua, u, v, cup, sup, bu)      ! Output: last three entries

         call get_east_boundary_field(nx, ny, u, ube)

         call get_east_boundary_field(nx, ny, v, vbe)

         call get_Ucbe_Vcbe(nx, ny, xe, ye, xke, yke, ube, vbe, Ucbe, Vcbe)

         call set_bcu(nx, ny, modvis, u, betan, gamman, au, bu) ! Output: last two entries

         call get_a9d_b_rescaling(nx, ny, au, bu)

         ! Calculating v coefficients and source
         call get_v_coefficients(nx, ny, coord, modvis, dt, rp, re, rn, Jp, &
            Je, Jn, xe, xk, alphae, betae, betan, gamman, vle, vln, vlp, roe, &
            ron, roa, Uce, Vcn, av) ! Output: last entry

         call get_v_source(nx, ny, coord, modvis, beta, dt, rp, re, rn &
            , xe, ye, xk, yk, xke, yke, xen, yen &
            , Jp, Je, Jn, roe, ron, roa, p, vle, vln &
            , Uce, Vcn, va, u, v, cvp, svp, bv) ! Output: last three entries

         call set_bcv(nx, ny, modvis, v, vin, betan, gamman, av, bv) ! Output: last two entries

         call get_a9d_b_rescaling(nx, ny, av, bv)

         ! Calculating SIMPLEC coefficients

         call get_internal_simplec_coefficients(nx, ny, re, rn, xe, ye, xk, &
            yk, au, av, due, dve, dun, dvn, de, dn) ! Last six are output

         call get_boundary_simplec_coefficients(nx, ny, de, dn) ! InOutput: de, dn

         g = 1.d0/(Rg*T)

         ! Calculating the coefficients of the linear system for pressure correction

         ! g, ro, Uce and Vcn used in this subroutine are those calculated in the previous iteraction
         ! de and dn must be calculated with the coef. of the linear sytem for u and v from which
         ! u*and v*are obtained.
         call get_p_coefficients(nx, ny, dt, rp, re, rn, Jp, Uce, Vcn, ro, g, &
            de, dn, ap) ! Output: last one

         ! Solving the linear system for u

         call lu2d9(au, dl9, du9, nxy, nx, ny) ! Matrix factorization

         call fb2d9(au, dl9, du9, ny, nx, nxy, u, tolu, nitm_u, bu) ! Solving factored linear system

         call norm_l1_9d(nx, ny, u, bu, au, norm)

         ! Solving the linear system for v

         call lu2d9(av, dl9, du9, nxy, nx, ny)  ! Matrix factorization

         call fb2d9(av, dl9, du9, ny, nx, nxy, v, tolu, nitm_u, bv) ! Solving factored linear system

         call norm_l1_9d(nx, ny, v, bv, av, norm)

         ! Calculating incorrect velocities at faces (Uce*and Vcn*)

         call get_velocities_at_faces(nx, ny, dt, rp, re, rn, xe, ye, xk, yk & ! Input
         , xen, yen, xke, yke, Jp, cup, cvp, sup, svp & ! Input
         , au, av, roa, p, u, v, uea, vea, una, vna & ! Input
         , ue, ve, un, vn, Uce, Vcn)                    ! Output

         call get_Uce_Vcn_at_boundary_faces(nx, ny, ye, u, Uce, Vcn) ! Output: Uce, Vcn

         ! Calculating source of the equation for pressure correction

         ! ro, Uce and Vcn are the incorrect ones (obtained with p*);
         ! rom, Ucem and Vcnm are those of the previous iteraction
         call get_p_source(nx, ny, beta, dt, rp, re, rn, Jp, ro, ro, roa, &
            Ucem, Uce, Vcnm, Vcn, bp) ! Output: last one

         ! Mass correction cycle

         do im = 1, imax

            call set_bcp(nx, ny, ap, bp) ! Output: last two entries

            call get_a5d_b_rescaling(nx, ny, ap, bp)

            call lu2d5(ap, dl5, du5, nxy, nx, ny) ! Matrix factorization

            call fb2d5(ap, dl5, du5, ny, nx, nxy, pl, tolp, nitm_p, bp) ! Solving factored linear system

         end do

         call norm_l1_5d(nx, ny, pl, bp, ap, norm)

         call get_pressure_density_correction_with_pl(nx, ny, pl, g, ro, p) ! InOutput: last two entries

         call get_p_extrapolation_to_fictitious(nx, ny, pin, p)

         call get_density_at_nodes(nx, ny, Rg, p, T, ro)

         call get_u_v_at_real_nodes_with_pl(nx, ny, xe, ye, xk, yk, rp, pl, &
            au, av, u, v) ! InOutput: last two entries

         call get_u_extrapolation_to_fictitious(nx, ny, modvis, betan, gamman, u) ! InOutput: last two entries

         call get_v_extrapolation_to_fictitious(nx, ny, modvis, betan, gamman, v) ! InOutput: last two entries

         call get_velocities_at_internal_faces_with_pl(nx, ny, xe, ye, xk, yk, &
            due, dve, dun, dvn, pl, ue, ve, un, vn, Uce, Vcn) ! Last six are InOutput

         call get_Uce_Vcn_at_boundary_faces_with_pl(nx, ny, pl, de, Uce) ! InOutput: Last entry

         call get_density_at_faces(nx, ny, beta, ro, Uce, Vcn, roe, ron) ! Output: last two entries

         ! Solving the linear system for temperature

         call get_T_coefficients_and_source(nx, ny, coord, modvis, beta, dt, &
            rp, re, rn, xe, ye, xk, yk, alphae, betae, betan, gamman, Jp, Je, &
            Jn, roe, ron, roa, p, pa, cp, vlp, ke, kn, Uce, Vcn, u, v, T, Ta, &
            at, bt) ! Output: last two entries

         call set_bcT(nx, ny, ccTw, Tin, Twall, T, betan, gamman, at, bt) ! Output: last two entries

         call lu2d9(at, dl9, du9, nxy, nx, ny) ! Matrix factorization

         call fb2d9(at, dl9, du9, ny, nx, nxy, T, tolu, nitm_u, bt) ! Solving factored linear system

         call norm_l1_9d(nx, ny, T, bt, at, norm)

         call get_density_at_nodes(nx, ny, Rg, p, T, ro) ! ro is output

         call get_density_at_faces(nx, ny, beta, ro, Uce, Vcn, roe, ron) ! roe and ron are output

         if (it < itdt) then
            ! Calculating convergence coefficients
            maxcu = get_maxc_9d(nx, ny, au)
            maxcv = get_maxc_9d(nx, ny, av)
            maxct = get_maxc_9d(nx, ny, at)
            maxcp = get_maxc_5d(nx, ny, ap)
            maxc = max(maxcu, maxcv, maxct, maxcp)
         end if

         ! Saving backup data
         if (mod(it, wbkp) == 0 .and. optm == 0) then

            ! Finishes cpu time
            call cpu_time(tcpu2)

            call save_backup(nx, ny, it, norm1, tcpuo + tcpu2 - tcpu1, p, pin, &
               pl, T, ro, roe, ron, u, v, ue, ve, un, vn, Uce, Vcn, pmod, Ns, &
               Yi) ! Output: none

         end if

         ! Printing results to listing file
         if (it == ito) then

            if (it == 1) norm1 = norm

            ! lid - listing file id
            write(lid,*)
            write(lid,"(A7, 4(A23))") "it", "norm/norm1", "fmi (kg/s)", &
               "fme (kg/s)", "Fd (N)"
            write(lid,*)

            write(*,*)
            write(*, "(A18, A7, A13, 3(A18))") "sim_id", "it", "norm/norm1", &
               "fmi (kg/s)", "fme (kg/s)", "Fd (N)"
            write(*,*)

            ! rid = residual file id
            write(rid, "('#', A5, 4(1X, A14))") "it", "DFd*", "norm/norm1", &
               "dt", "maxc"

         end if

         ! Printing results to listing file
         if (mod(it, wlf) == 0) then

            call get_mass_flow_rate_and_thrust(nx, ny, re, roe, u, Uce, fmi, &
               fme, Fd) ! Output: last three entries

            dFd = dabs(Fd/Fd1D - Fda_a)

            write(lid,"(I7, 4(ES23.15))") it, norm/norm1, fmi, fme, Fd

            if (len(trim(adjustl(sim_id))) <= 18) then
               write(*,73) trim(adjustl(sim_id)), it, norm/norm1, fmi, fme, Fd
               73 format(A18, I7, ES13.5, 3(ES18.10))
            else
               write(*,74) trim(adjustl(sim_id)), it, norm/norm1, fmi, fme, Fd
               74 format(A, I7, ES13.5, 3(ES18.10))
            end if

            write(rid, "(I6, 4(1X, ES14.7))") it, dFd, norm/norm1, dt, maxc

            Fda_a = Fd/Fd1D

            call flush(lid)
            call flush(rid)
         end if

         if (it < itdt) call get_time_step_increment(maxc, maxcc, dt)

         ! Parada com base no valor 'não definido' da variável norm/norm1
         if (norm/norm1 .NE. norm/norm1) then

            write(*,51)
            write(lid,51)
            write(rid,51)
            flush(lid)
            flush(rid)
            close(lid)
            close(rid)

            if (optm /= 0) then

               call depp_save_fitness(0.d0, 1, "ERROR: divergence")

            end if

            stop

         endif

         if (optm /= 0 .and. it == itmax) then

            call depp_save_fitness(0.d0, 0, "ERROR: ITMAX reached")

            stop

         end if

      end do

   end subroutine get_one_species_constant

   !============================================================================

   subroutine get_one_species_variable

      implicit none
      integer :: im
      integer :: cont = 0
      real*8 :: dFd

      51 format(" ==========================================================", &
         "==================== ", /, " ============ ERROR: Numero nao ", &
         "definido - Tente outros parametros ============ ", /, &
         " =================================================================", &
         "============= ", /)

      do it = ito, itmax

         ! Parada com base na tolerância para o dFd*
         if (cont > ittol) then
            exit
         else
            if (dFd < tolerance) then
               cont = cont + 1
            else
               cont = 0
            end if
         endif

         norm = 0.d0

         ! Calculating mean variables used in the eddy viscosity calculation
         if (modvis == 1 .and. modtur == 1) then
            ro_av = 0.5d0*(ro + roa)
            ron_av = 0.5d0*(ron + rona)
            u_av = 0.5d0*(u + ua)
            v_av = 0.5d0*(v + va)
         end if

         ! Updating variables for the new time step
         ua = u
         va = v
         uea = ue
         vea = ve
         una = un
         vna = vn
         pa = p
         Ta = T
         roa = ro
         rona = ron
         Ucem = Uce
         Vcnm = Vcn
         pina = pin

         ! Updating conditions at the nozzle entrance
         call get_uin_vin_pin_Tin_Mw(nx, ny, gamma, Rg, po, T0, u & ! Input
         , uin, vin, pin, Tin, Mw)        ! Output

         call get_plin_and_p_fictitious(nx, ny, pina, pin, plin, p) ! Output: last two entries

         ! Calculation of the gas properties
         call set_variable_cp_and_gamma(folder_input, nx, ny, specie, Rg, T, &
            cp, gcp) ! Output: last two entries

         if (modvis == 1) then

            call set_variable_laminar_viscosity_at_nodes(folder_input, nx, ny, &
               specie, T, vlp) ! Output: last entry

            call get_laminar_viscosity_at_faces(nx, ny, vlp, vle, vln) ! Output: last two entries

            if (modtur == 1) then

               ! Calculates the dimensional eddy viscosity at center of real volumes
               call get_eddy_viscosity_blomax(nx, ny, x, y, xp, yp, xk, yk, &
                  Jp, Jn, vln, ro_av, ron_av, u_av, v_av, vtp) ! Last one is output

               ! Extrapolation of the eddy viscosity to the fictitious volumes
               call get_eddy_viscosity_fictitious(nx, ny, vtp)

               vlp = vlp + vtp

               ! CAUTION: In order to apply correctly this function, laminar viscosity vlp must be known
               ! in real and fictitious volumes (except CV of corners SW, SE, NW and NE)
               call get_laminar_viscosity_at_faces(nx, ny, vlp, vle, vln) ! Output: last two entries

            end if

            call set_variable_thermal_conductivity_at_nodes(folder_input, nx, &
               ny, specie, T, kp) ! Output: last entry

            if (modtur == 1) then
               kp = kp + vtp*cp/0.9d0
            end if

            call get_thermal_conductivity_at_faces(nx, ny, kp, ke, kn) ! Output: last two entries

         end if

         ! Calculating u coefficients and source
         call get_u_coefficients(nx, ny, modvis, dt, rp, re, rn, Jp, Je, Jn &
            , ye, yk, alphae, betae, betan, gamman &
            , vle, vln, roe, ron, roa, Uce, Vcn, au)    ! Output: last entry

         call get_u_source(nx, ny, modvis, beta, dt, rp, re, rn&
            , xe, ye, xk, yk, xke, yke, xen, yen &
            , Jp, Je, Jn, roe, ron, roa, p, vle, vln &
            , Uce, Vcn, ua, u, v, cup, sup, bu)      ! Output: last three entries

         call get_east_boundary_field(nx, ny, u, ube)

         call get_east_boundary_field(nx, ny, v, vbe)

         call get_Ucbe_Vcbe(nx, ny, xe, ye, xke, yke, ube, vbe, Ucbe, Vcbe)

         call set_bcu(nx, ny, modvis, u, betan, gamman, au, bu) ! Output: last two entries

         call get_a9d_b_rescaling(nx, ny, au, bu)

         ! Calculating v coefficients and source
         call get_v_coefficients(nx, ny, coord, modvis, dt, rp, re, rn, Jp, &
            Je, Jn, xe, xk, alphae, betae, betan, gamman, vle, vln, vlp, roe, &
            ron, roa, Uce, Vcn, av) ! Output: last entry

         call get_v_source(nx, ny, coord, modvis, beta, dt, rp, re, rn &
            , xe, ye, xk, yk, xke, yke, xen, yen &
            , Jp, Je, Jn, roe, ron, roa, p, vle, vln &
            , Uce, Vcn, va, u, v, cvp, svp, bv) ! Output: last three entries

         call set_bcv(nx, ny, modvis, v, vin, betan, gamman, av, bv) ! Output: last two entries

         call get_a9d_b_rescaling(nx, ny, av, bv)

         ! Calculating SIMPLEC coefficients

         call get_internal_simplec_coefficients(nx, ny, re, rn, xe, ye, xk, &
            yk, au, av, due, dve, dun, dvn, de, dn) ! Last six are Output

         call get_boundary_simplec_coefficients(nx, ny, de, dn) ! InOutput: de, dn

         g = 1.d0/(Rg*T)

         ! Calculating the coefficients of the linear system for pressure correction

         ! g, ro, Uce and Vcn used in this subroutine are those calculated in the previous iteraction
         ! de and dn must be calculated with the coef. of the linear sytem for u and v from which
         ! u*and v*are obtained.
         call get_p_coefficients(nx, ny, dt, rp, re, rn, Jp, Uce, Vcn, ro, g, &
            de, dn, ap) ! Output: last one

         ! Solving the linear system for u

         call lu2d9(au, dl9, du9, nxy, nx, ny) ! Matrix factorization

         call fb2d9(au, dl9, du9, ny, nx, nxy, u, tolu, nitm_u, bu) ! Solving factored linear system

         call norm_l1_9d(nx, ny, u, bu, au, norm)

         ! Solving the linear system for v

         call lu2d9(av, dl9, du9, nxy, nx, ny)  ! Matrix factorization

         call fb2d9(av, dl9, du9, ny, nx, nxy, v, tolu, nitm_u, bv) ! Solving factored linear system

         call norm_l1_9d(nx, ny, v, bv, av, norm)

         ! Calculating incorrect velocities at faces (Uce*and Vcn*)

         call get_velocities_at_faces(nx, ny, dt, rp, re, rn, xe, ye, xk, yk & ! Input
         , xen, yen, xke, yke, Jp, cup, cvp, sup, svp & ! Input
         , au, av, roa, p, u, v, uea, vea, una, vna & ! Input
         , ue, ve, un, vn, Uce, Vcn)                    ! Output

         call get_Uce_Vcn_at_boundary_faces(nx, ny, ye, u, Uce, Vcn) ! Output: Uce, Vcn

         ! Calculating source of the equation for pressure correction

         ! ro, Uce and Vcn are the incorrect ones (obtained with p*);
         ! rom, Ucem and Vcnm are those of the previous iteraction
         call get_p_source(nx, ny, beta, dt, rp, re, rn, Jp, ro, ro, roa, &
            Ucem, Uce, Vcnm, Vcn, bp) ! Output: last one

         ! Mass correction cycle

         do im = 1, imax

            call set_bcp(nx, ny, ap, bp) ! Output: last two entries

            call get_a5d_b_rescaling(nx, ny, ap, bp)

            call lu2d5(ap, dl5, du5, nxy, nx, ny) ! Matrix factorization

            call fb2d5(ap, dl5, du5, ny, nx, nxy, pl, tolp, nitm_p, bp) ! Solving factored linear system

         end do

         call norm_l1_5d(nx, ny, pl, bp, ap, norm)

         call get_pressure_density_correction_with_pl(nx, ny, pl, g, ro, p) ! InOutput: last two entries

         call get_p_extrapolation_to_fictitious(nx, ny, pin, p)

         call get_density_at_nodes(nx, ny, Rg, p, T, ro)

         call get_u_v_at_real_nodes_with_pl(nx, ny, xe, ye, xk, yk, rp, pl, &
            au, av, u, v) ! InOutput: last two entries

         call get_u_extrapolation_to_fictitious(nx, ny, modvis, betan, gamman, u) ! InOutput: last two entries

         call get_v_extrapolation_to_fictitious(nx, ny, modvis, betan, gamman, v) ! InOutput: last two entries

         call get_velocities_at_internal_faces_with_pl(nx, ny, xe, ye, xk, yk, &
            due, dve, dun, dvn, pl, ue, ve, un, vn, Uce, Vcn) ! Last six are InOutput

         call get_Uce_Vcn_at_boundary_faces_with_pl(nx, ny, pl, de, Uce) ! InOutput: Last entry

         call get_density_at_faces(nx, ny, beta, ro, Uce, Vcn, roe, ron) ! Output: last two entries

         ! Solving the linear system for temperature

         call get_T_coefficients_and_source(nx, ny, coord, modvis, beta, dt, &
            rp, re, rn, xe, ye, xk, yk, alphae, betae, betan, gamman, Jp, Je, &
            Jn, roe, ron, roa, p, pa, cp, vlp, ke, kn, Uce, Vcn, u, v, T, Ta, &
            at, bt) ! Output: last two entries

         call set_bcT(nx, ny, ccTw, Tin, Twall, T, betan, gamman, at, bt) ! Output: last two entries

         call lu2d9(at, dl9, du9, nxy, nx, ny) ! Matrix factorization

         call fb2d9(at, dl9, du9, ny, nx, nxy, T, tolu, nitm_u, bt) ! Solving factored linear system

         call norm_l1_9d(nx, ny, T, bt, at, norm)

         call get_density_at_nodes(nx, ny, Rg, p, T, ro) ! ro is output

         call get_density_at_faces(nx, ny, beta, ro, Uce, Vcn, roe, ron) ! roe and ron are output

         if (it < itdt) then
            ! Calculating convergence coefficients
            maxcu = get_maxc_9d(nx, ny, au)
            maxcv = get_maxc_9d(nx, ny, av)
            maxct = get_maxc_9d(nx, ny, at)
            maxcp = get_maxc_5d(nx, ny, ap)
            maxc = max(maxcu, maxcv, maxct, maxcp)
         end if

         ! Saving backup data
         if (mod(it, wbkp) == 0 .and. optm == 0) then

            ! Finishes cpu time
            call cpu_time(tcpu2)

            call save_backup(nx, ny, it, norm1, tcpuo + tcpu2 - tcpu1, p, pin, &
               pl, T, ro, roe, ron, u, v, ue, ve, un, vn, Uce, Vcn, pmod, &
               Ns, Yi) ! Output: none

         end if

         ! Printing results to listing file
         if (it == ito) then

            if (it == 1) norm1 = norm

            ! lid - listing file id
            write(lid,*)
            write(lid,"(A7, 4(A23))") "it", "norm/norm1", "fmi (kg/s)", &
               "fme (kg/s)", "Fd (N)"
            write(lid,*)

            write(*,*)
            write(*, "(A18, A7, A13, 3(A18))") "sim_id", "it", "norm/norm1", &
               "fmi (kg/s)", "fme (kg/s)", "Fd (N)"
            write(*,*)

            ! rid = residual file id
            write(rid, "('#', A5, 4(1X, A14))") "it", "DFd*", "norm/norm1", &
               "dt", "maxc"

         end if

         ! Printing results to listing file
         if (mod(it, wlf) == 0) then

            call get_mass_flow_rate_and_thrust(nx, ny, re, roe, u, Uce, fmi, &
               fme, Fd) ! Output: last three entries

            dFd = dabs(Fd/Fd1D - Fda_a)

            write(lid,"(I7, 4(ES23.15))") it, norm/norm1, fmi, fme, Fd

            if (len(trim(adjustl(sim_id))) <= 18) then
               write(*,73) trim(adjustl(sim_id)), it, norm/norm1, fmi, fme, Fd
               73 format(A18, I7, ES13.5, 3(ES18.10))
            else
               write(*,74) trim(adjustl(sim_id)), it, norm/norm1, fmi, fme, Fd
               74 format(A, I7, ES13.5, 3(ES18.10))
            end if

            write(rid, "(I6, 4(1X, ES14.7))") it, dFd, norm/norm1, dt, maxc

            Fda_a = Fd/Fd1D

            call flush(lid)
            call flush(rid)
         end if

         if (it < itdt) call get_time_step_increment(maxc, maxcc, dt)

         ! Parada com base no valor 'não definido' da variável norm/norm1
         if (norm/norm1 .NE. norm/norm1) then

            write(*,51)
            write(lid,51)
            write(rid,51)
            flush(lid)
            flush(rid)
            close(lid)
            close(rid)

            if (optm /= 0) then

               call depp_save_fitness(0.d0, 1, "ERROR: divergence")

            end if

            stop

         endif

         if (optm /= 0 .and. it == itmax) then

            call depp_save_fitness(0.d0, 0, "ERROR: ITMAX reached")

            stop

         end if

      end do

   end subroutine get_one_species_variable

   !============================================================================

   subroutine get_frozen_flow

      implicit none
      integer :: i, j, np, im
      integer :: unit_au = 82 ! identifier for the auxiliary fields file
      integer :: cont = 0
      real*8 :: dFd

      real(8), allocatable, dimension(:) :: temp_aid  ! auxiliar vector for temperature at wall

      51 format(" ==========================================================", &
         "==================== ", /, " ============ ERROR: Numero nao ", &
         "definido - Tente outros parametros ============ ", /, &
         " =================================================================", &
         "============= ", /)

      if (cp_type == 2) then
         call GIBBS_EQUILIBRIO_dados_recebe(unit_au, 0, cmod, itimax, itemax, &
            tol_ej, tol_n, T0, po/101325.0d0, OF)
         call GIBBS_EQUILIBRIO_inicio1(folder_input)
         call GIBBS_EQUILIBRIO_calculos1
         call GIBBS_MACH1D_envia2 (Ns)
         allocate(Yi(nx*ny, Ns), species(Ns))
         call GIBBS_MACH1D_envia3 (species)
         call GIBBS_MACH1D_envia6 (cp_cte, gamma, Rg, Yi(1, :))
         do i = 2, nx*ny
            Yi(i, :) = Yi(1, :)
         end do
         call frozen_cp_gamma (nx, ny, T, Yi, cmod, itimax, itemax, tol_ej, &
            tol_n, po, OF, cp, gcp)
      else
         allocate(Yi(nx*ny,1), species(1))
         Yi = 0.0d0
         Ns = 1
         species = 'none'
         cp = cp_cte
         gcp = gamma
      end if
      Rgp = Rg

      ! Properties - initialization
      call GIBBS_mi_kc_inicio(folder_input)
      call frozen_mi_k(nx, ny, T, v_type, visc_cte, k_type, k_cte, vlp, kp)
      call frozen_Mach(Rgp, gcp, u, T, M)

      ! Temporal cycle
      do it = ito, itmax

         ! Parada com base na tolerância para o dFd*
         if (cont > ittol) then
            exit
         else
            if (dFd < tolerance) then
               cont = cont + 1
            else
               cont = 0
            end if
         endif

         norm = 0.d0

         ! Calculating mean variables used in the eddy viscosity calculation
         if (modvis == 1 .and. modtur == 1) then
            ro_av = 0.5d0*(ro + roa)
            ron_av = 0.5d0*(ron + rona)
            u_av = 0.5d0*(u + ua)
            v_av = 0.5d0*(v + va)
         end if

         ! Updating variables for the new time step
         ua = u
         va = v
         uea = ue
         vea = ve
         una = un
         vna = vn
         pa = p
         Ta = T
         roa = ro
         rona = ron
         Ucem = Uce
         Vcnm = Vcn
         pina = pin

         ! Updating conditions at the nozzle entrance
         call get_uin_vin_pin_Tin_Mw(nx, ny, gamma, Rg, po, T0, u & ! Input
         , uin, vin, pin, Tin, Mw)        ! Output

         call get_plin_and_p_fictitious(nx, ny, pina, pin, plin, p) ! Output: last two entries

         ! Evaluation of chemical properties
         if ((it >= it_ch) .and. ((it == 1) .or. (it == itmax) .or. &
               (mod(it, w_ch) == 0))) then

            ! Numerical scheme to avoid the computation of negative temperatures (Part 1: prescribed wall temperature)
            if (ccTw == 1) then
               if (.not. allocated(temp_aid)) allocate(temp_aid(Nx))
               j = ny
               do i = 1, nx
                  np = (j-1)*nx + i
                  temp_aid(i) = T(np)
                  T(np) = Twall(i)
               end do
            end if

            if (cp_type == 2) then
               call GIBBS_EQUILIBRIO_dados_recebe (unit_au, 0, cmod, itimax, &
                  itemax, tol_ej, tol_n, T0, po/101325.0d0, OF)
               call GIBBS_EQUILIBRIO_calculos1
               call GIBBS_MACH1D_envia2 (Ns)
               call GIBBS_MACH1D_envia3 (species)
               do j = 1, ny
                  call GIBBS_MACH1D_envia6 (cp_cte, gamma, Rg, Yi(j, :))
               end do
               do np = 2, nxy
                  Yi(np,:) = Yi(1, :)
               end do

               call frozen_cp_gamma(nx, ny, T, Yi, cmod, itimax, itemax, &
                  tol_ej, tol_n, po, OF, cp, gcp) ! Output: last two entries
            else
               allocate(Yi(nxy,1), species(1))
               Yi = 0.0d0
               Ns = 1
               species = 'none'
               cp = cp_cte
               gcp = gamma
            end if

            call frozen_mi_k(nx, ny, T, v_type, visc_cte, k_type, k_cte, vlp, &
               kp)  ! Output: last two entries

            ! Numerical scheme to avoid the computation of negative temperatures (Part 2: return)
            if (ccTw == 1) then
               j = ny
               do i = 1, nx
                  np = (j-1)*nx + i
                  T(np) = temp_aid(i)
               end do
            end if

         end if

         if (modvis == 1) then

            ! Calculation of the gas properties
            call get_laminar_viscosity_at_faces(nx, ny, vlp, vle, vln) ! Output: last two entries

            if (modtur == 1) then

               ! Calculates the dimensional eddy viscosity at center of real volumes
               call get_eddy_viscosity_blomax(nx, ny, x, y, xp, yp, xk, yk, &
                  Jp, Jn, vln, ro_av, ron_av, u_av, v_av, vtp) ! Last one is output

               ! Extrapolation of the eddy viscosity to the fictitious volumes
               call get_eddy_viscosity_fictitious(nx, ny, vtp)

               vlp = vlp + vtp

               ! CAUTION: In order to apply correctly this function, laminar viscosity vlp must be known
               ! in real and fictitious volumes (except CV of corners SW, SE, NW and NE)
               call get_laminar_viscosity_at_faces(nx, ny, vlp, vle, vln) ! Output: last two entries

            end if

            if (modtur == 1) then
               kp = kp + vtp*cp/0.9d0
            end if

            call get_thermal_conductivity_at_faces(nx, ny, kp, ke, kn) ! Output: last two entries

         end if

         ! Calculating u coefficients and source
         call get_u_coefficients(nx, ny, modvis, dt, rp, re, rn, Jp, Je, Jn &
            , ye, yk, alphae, betae, betan, gamman &
            , vle, vln, roe, ron, roa, Uce, Vcn, au)    ! Output: last entry

         call get_u_source(nx, ny, modvis, beta, dt, rp, re, rn&
            , xe, ye, xk, yk, xke, yke, xen, yen &
            , Jp, Je, Jn, roe, ron, roa, p, vle, vln &
            , Uce, Vcn, ua, u, v, cup, sup, bu)      ! Output: last three entries

         call get_east_boundary_field(nx, ny, u, ube)

         call get_east_boundary_field(nx, ny, v, vbe)

         call get_Ucbe_Vcbe(nx, ny, xe, ye, xke, yke, ube, vbe, Ucbe, Vcbe)

         call set_bcu(nx, ny, modvis, u, betan, gamman, au, bu) ! Output: last two entries

         call get_a9d_b_rescaling(nx, ny, au, bu)

         ! Calculating v coefficients and source
         call get_v_coefficients(nx, ny, coord, modvis, dt, rp, re, rn, Jp, &
            Je, Jn, xe, xk, alphae, betae, betan, gamman, vle, vln, vlp, roe, &
            ron, roa, Uce, Vcn, av) ! Output: last entry

         call get_v_source(nx, ny, coord, modvis, beta, dt, rp, re, rn &
            , xe, ye, xk, yk, xke, yke, xen, yen &
            , Jp, Je, Jn, roe, ron, roa, p, vle, vln &
            , Uce, Vcn, va, u, v, cvp, svp, bv) ! Output: last three entries

         call set_bcv(nx, ny, modvis, v, vin, betan, gamman, av, bv) ! Output: last two entries

         call get_a9d_b_rescaling(nx, ny, av, bv)

         ! Calculating SIMPLEC coefficients

         call get_internal_simplec_coefficients(nx, ny, re, rn, xe, ye, xk, &
            yk, au, av, due, dve, dun, dvn, de, dn) ! Last six are Output

         call get_boundary_simplec_coefficients(nx, ny, de, dn) ! InOutput: de, dn

         g = 1.d0/(Rg*T)

         ! Calculating the coefficients of the linear system for pressure correction

         ! g, ro, Uce and Vcn used in this subroutine are those calculated in the previous iteraction
         ! de and dn must be calculated with the coef. of the linear sytem for u and v from which
         ! u*and v*are obtained.
         call get_p_coefficients(nx, ny, dt, rp, re, rn, Jp, Uce, Vcn, ro, g, &
            de, dn, ap) ! Output: last one

         ! Solving the linear system for u

         call lu2d9(au, dl9, du9, nxy, nx, ny) ! Matrix factorization

         call fb2d9(au, dl9, du9, ny, nx, nxy, u, tolu, nitm_u, bu) ! Solving factored linear system

         call norm_l1_9d(nx, ny, u, bu, au, norm)

         ! Solving the linear system for v

         call lu2d9(av, dl9, du9, nxy, nx, ny)  ! Matrix factorization

         call fb2d9(av, dl9, du9, ny, nx, nxy, v, tolu, nitm_u, bv) ! Solving factored linear system

         call norm_l1_9d(nx, ny, v, bv, av, norm)

         ! Calculating incorrect velocities at faces (Uce*and Vcn*)

         call get_velocities_at_faces(nx, ny, dt, rp, re, rn, xe, ye, xk, yk & ! Input
         , xen, yen, xke, yke, Jp, cup, cvp, sup, svp & ! Input
         , au, av, roa, p, u, v, uea, vea, una, vna & ! Input
         , ue, ve, un, vn, Uce, Vcn)                    ! Output

         call get_Uce_Vcn_at_boundary_faces(nx, ny, ye, u, Uce, Vcn) ! Output: Uce, Vcn

         ! Calculating source of the equation for pressure correction

         ! ro, Uce and Vcn are the incorrect ones (obtained with p*);
         ! rom, Ucem and Vcnm are those of the previous iteraction
         call get_p_source(nx, ny, beta, dt, rp, re, rn, Jp, ro, ro, roa, &
            Ucem, Uce, Vcnm, Vcn, bp) ! Output: last one

         ! Mass correction cycle

         do im = 1, imax

            call set_bcp(nx, ny, ap, bp) ! Output: last two entries

            call get_a5d_b_rescaling(nx, ny, ap, bp)

            call lu2d5(ap, dl5, du5, nxy, nx, ny) ! Matrix factorization

            call fb2d5(ap, dl5, du5, ny, nx, nxy, pl, tolp, nitm_p, bp) ! Solving factored linear system

         end do

         call norm_l1_5d(nx, ny, pl, bp, ap, norm)

         call get_pressure_density_correction_with_pl(nx, ny, pl, g, ro, p) ! InOutput: last two entries

         call get_p_extrapolation_to_fictitious(nx, ny, pin, p)

         call get_density_at_nodes(nx, ny, Rg, p, T, ro)

         call get_u_v_at_real_nodes_with_pl(nx, ny, xe, ye, xk, yk, rp, pl, &
            au, av, u, v) ! InOutput: last two entries

         call get_u_extrapolation_to_fictitious(nx, ny, modvis, betan, gamman, u) ! InOutput: last two entries

         call get_v_extrapolation_to_fictitious(nx, ny, modvis, betan, gamman, v) ! InOutput: last two entries

         call get_velocities_at_internal_faces_with_pl(nx, ny, xe, ye, xk, yk, &
            due, dve, dun, dvn, pl, ue, ve, un, vn, Uce, Vcn) ! Last six are InOutput

         call get_Uce_Vcn_at_boundary_faces_with_pl(nx, ny, pl, de, Uce) ! InOutput: Last entry

         call get_density_at_faces(nx, ny, beta, ro, Uce, Vcn, roe, ron) ! Output: last two entries

         ! Solving the linear system for temperature

         call get_T_coefficients_and_source(nx, ny, coord, modvis, beta, dt, &
            rp, re, rn, xe, ye, xk, yk, alphae, betae, betan, gamman, Jp, Je, &
            Jn, roe, ron, roa, p, pa, cp, vlp, ke, kn, Uce, Vcn, u, v, T, Ta, &
            at, bt) ! Output: last two entries

         call set_bcT(nx, ny, ccTw, Tin, Twall, T, betan, gamman, at, bt) ! Output: last two entries

         call lu2d9(at, dl9, du9, nxy, nx, ny) ! Matrix factorization

         call fb2d9(at, dl9, du9, ny, nx, nxy, T, tolu, nitm_u, bt) ! Solving factored linear system

         call norm_l1_9d(nx, ny, T, bt, at, norm)

         call get_density_at_nodes(nx, ny, Rg, p, T, ro) ! ro is output

         call get_density_at_faces(nx, ny, beta, ro, Uce, Vcn, roe, ron) ! roe and ron are output

         if (it < itdt) then
            ! Calculating convergence coefficients
            maxcu = get_maxc_9d(nx, ny, au)
            maxcv = get_maxc_9d(nx, ny, av)
            maxct = get_maxc_9d(nx, ny, at)
            maxcp = get_maxc_5d(nx, ny, ap)
            maxc = max(maxcu, maxcv, maxct, maxcp)
         end if

         ! Saving backup data
         if (mod(it, wbkp) == 0 .and. optm == 0) then

            ! Finishes cpu time
            call cpu_time(tcpu2)

            call save_backup(nx, ny, it, norm1, tcpuo + tcpu2 - tcpu1, p, pin, &
               pl, T, ro, roe, ron, u, v, ue, ve, un, vn, Uce, Vcn, pmod, &
               Ns, Yi) ! Output: none

         end if

         ! Printing results to listing file
         if (it == ito) then

            if (it == 1) norm1 = norm

            ! lid - listing file id
            write(lid,*)
            write(lid,"(A7, 4(A23))") "it", "norm/norm1", "fmi (kg/s)", &
               "fme (kg/s)", "Fd (N)"
            write(lid,*)

            write(*,*)
            write(*, "(A18, A7, A13, 3(A18))") "sim_id", "it", "norm/norm1", &
               "fmi (kg/s)", "fme (kg/s)", "Fd (N)"
            write(*,*)

            ! rid = residual file id
            write(rid, "('#', A5, 4(1X, A14))") "it", "DFd*", "norm/norm1", &
               "dt", "maxc"

         end if

         ! Printing results to listing file
         if (mod(it, wlf) == 0) then

            call get_mass_flow_rate_and_thrust(nx, ny, re, roe, u, Uce, fmi, &
               fme, Fd) ! Output: last three entries

            dFd = dabs(Fd/Fd1D - Fda_a)

            write(lid,"(I7, 4(ES23.15))") it, norm/norm1, fmi, fme, Fd

            if (len(trim(adjustl(sim_id))) <= 18) then
               write(*,73) trim(adjustl(sim_id)), it, norm/norm1, fmi, fme, Fd
               73 format(A18, I7, ES13.5, 3(ES18.10))
            else
               write(*,74) trim(adjustl(sim_id)), it, norm/norm1, fmi, fme, Fd
               74 format(A, I7, ES13.5, 3(ES18.10))
            end if

            write(rid, "(I6, 4(1X, ES14.7))") it, dFd, norm/norm1, dt, maxc

            Fda_a = Fd/Fd1D

            call flush(lid)
            call flush(rid)
         end if

         if (it < itdt) call get_time_step_increment(maxc, maxcc, dt)

         ! Parada com base no valor 'não definido' da variável norm/norm1
         if (norm/norm1 .NE. norm/norm1) then

            write(*,51)
            write(lid,51)
            write(rid,51)
            flush(lid)
            flush(rid)
            close(lid)
            close(rid)

            if (optm /= 0) then

               call depp_save_fitness(0.d0, 1, "ERROR: divergence")

            end if

            stop

         endif

         if (optm /= 0 .and. it == itmax) then

            call depp_save_fitness(0.d0, 0, "ERROR: ITMAX reached")

            stop

         end if

      end do

   end subroutine get_frozen_flow

   !============================================================================

   subroutine get_local_equilibrium

      integer :: i, j, ii           ! counters and indeces
      integer :: im                 ! counters - iterations
      integer :: np                 ! counter: control volume
      integer :: npe, npee          ! neighbour: east, east-east
      integer :: npw, npww          ! neighbour: west, west-west
      integer :: unit_au = 82       ! identifier for the auxiliary fields file
      integer :: cont = 0
      real*8 :: dFd
      real(8) :: dxf, dxr           ! auxiliary variables for inlet nozzle properties
      real(8) :: x_e, x_w           ! auxiliary variables for distances

      real(8), allocatable, dimension(:) :: temp_aid  ! auxiliar vector for temperature at wall

      51 format(" ==========================================================", &
         "==================== ", /, " ============ ERROR: Numero nao ", &
         "definido - Tente outros parametros ============ ", /, &
         " =================================================================", &
         "============= ", /)

      ! Initial chemical/physical properties
      if (cp_type == 2) then
         call GIBBS_EQUILIBRIO_dados_recebe (unit_au, 0, cmod, itimax, itemax, &
            tol_ej, tol_n, T0, po/101325.0d0, OF)
         call GIBBS_EQUILIBRIO_inicio1(folder_input)
         call GIBBS_EQUILIBRIO_calculos1
         call GIBBS_MACH1D_envia2 (Ns)
         call GIBBS_MACH1D_envia9 (Nr)
         allocate(Yi(nxy, Ns), Yia(nxy, Ns), species(Ns), hi(nxy, Ns), &
            gcpf(nxy), entin(ny), ent(nxy), ej(nxy, Nr))
         call GIBBS_MACH1D_inicio1
         call GIBBS_hi_massa_calculo
         call GIBBS_MACH1D_envia3 (species)
         call GIBBS_MACH1D_envia4 (cp_cte, gamma, Rg, Yi(1, :), hi(1, :))
         call eq_frozen_cp_gamma(nx, ny, T, p, cmod, Nr, it, itmax, itimax, &
            itemax, tol_ej, tol_n, OF, Yi, hi, Rgp, cp, gcp, ej)
         call GIBBS_mi_kc_inicio(folder_input)
         gcpf = gcp
         call eq_effective_gamma(nx, ny, T, p, cmod, itimax, itemax, tol_ej, &
            tol_n, OF, dpY, Rgp, gcpf, gcp)
         Yia = 0.0d0
      else
         allocate(Yi(nxy,1), Yia(nxy, Ns), species(1), hi(nxy,1))
         Yi = 0.0d0
         Yia = 0.0d0
         Ns = 1
         species = 'none'
         cp = cp_cte
         gcp = gamma
         Rgp = Rg
      end if

      ! Evalution of Yi in ghost cells according to inlet conditions
      i = 1
      do j = 2, ny-1
         np = (j-1)*nx + i
         npe = np + 1
         npee = np + 2
         x_e = x(np)
         dxf = 2*(xp(npe) - x_e)
         dxr = xp(npee) - xp(npe)
         do ii = 1, Ns
            Yi(np,ii) = Yi(npe,ii) + (dxf/dxr)*(Yi(npe,ii) - Yi(npee, ii))
         end do
      end do

      ! Evaluation of Yi in ghost cells according to outlet conditions
      i = nx
      do j = 2, ny-1
         np = (j-1)*nx + i
         npw = np - 1
         npww = np - 2
         x_w = x(npw)
         dxf = 2*(x_w - xp(npw))
         dxr = xp(npw) - xp(npww)
         do ii = 1, Ns
            Yi(np,ii) = Yi(npw,ii) + (dxf/dxr)*(Yi(npw,ii) - Yi(npww, ii))
         end do
      end do

      ! Mach number initialization
      call eq_Mach (Rgp, gcp, u, T, M)

      ! Iterative process
      do it = ito, itmax

         ! Parada com base na tolerância para o dFd*
         if (cont > ittol) then
            exit
         else
            if (dFd < tolerance) then
               cont = cont + 1
            else
               cont = 0
            end if
         endif

         norm = 0.d0

         ! Calculating mean variables used in the eddy viscosity calculation
         if (modvis == 1 .and. modtur == 1) then
            ro_av = 0.5d0*(ro + roa)
            ron_av = 0.5d0*(ron + rona)
            u_av = 0.5d0*(u + ua)
            v_av = 0.5d0*(v + va)
         end if

         ! Updating variables for the new time step
         ua = u
         va = v
         uea = ue
         vea = ve
         una = un
         vna = vn
         pa = p
         Ta = T
         roa = ro
         rona = ron
         Ucem = Uce
         Vcnm = Vcn
         pina = pin

         if (it > 1) Yia = Yi

         ! Updating conditions at the nozzle entrance
         call get_uin_vin_pin_Tin_Mw(nx, ny, gamma, Rg, po, T0, u & ! Input
         , uin, vin, pin, Tin, Mw)        ! Output

         call get_plin_and_p_fictitious(nx, ny, pina, pin, plin, p) ! Output: last two entries

         ! Evaluation of chemical properties
         if ((it >= it_ch) .and. ((it == 1) .or. (it == itmax) .or. &
               (mod(it, w_ch) == 0))) then

            ! Numerical scheme to avoid the computation of negative temperatures (Part 1: prescribed wall temperature)
            if (ccTw == 1) then
               if (.not. allocated(temp_aid)) allocate(temp_aid(Nx))
               j = ny
               do i = 1, nx
                  np = (j-1)*nx + i
                  temp_aid(i) = T(np)
                  T(np) = Twall(i)
               end do
            end if

            if (cp_type == 2) call eq_frozen_cp_gamma(nx, ny, T, p, cmod, Nr, &
               it, itmax, itimax, itemax, tol_ej, tol_n, OF, Yi, hi, Rgp, &
               cp, gcp, ej) ! Output: the last three entries

            call eq_frozen_mi_k(nx, ny, T, v_type, visc_cte, k_type, k_cte, &
               vlp, kp)  ! Output: last two entries

            ! Numerical scheme to avoid the computation of negative temperatures (Part 2: return)
            if (ccTw == 1) then
               j = ny
               do i = 1, nx
                  np = (j-1)*nx + i
                  T(np) = temp_aid(i)
               end do
            end if

            if (modvis == 1) then

               ! Calculation of the gas properties
               call get_laminar_viscosity_at_faces(nx, ny, vlp, vle, vln) ! Output: last two entries

               if (modtur == 1) then

                  ! Calculates the dimensional eddy viscosity at center of real volumes
                  call get_eddy_viscosity_blomax(nx, ny, x, y, xp, yp, xk, yk, &
                     Jp, Jn, vln, ro_av, ron_av, u_av, v_av, vtp) ! Last one is output

                  ! Extrapolation of the eddy viscosity to the fictitious volumes
                  call get_eddy_viscosity_fictitious(nx, ny, vtp)

                  vlp = vlp + vtp

                  ! CAUTION: In order to apply correctly this function, laminar viscosity vlp must be known
                  ! in real and fictitious volumes (except CV of corners SW, SE, NW and NE)
                  call get_laminar_viscosity_at_faces(nx, ny, vlp, vle, vln) ! Output: last two entries

               end if

               if (modtur == 1) then
                  kp = kp + vtp*cp/0.9d0
               end if

               call get_thermal_conductivity_at_faces(nx, ny, kp, ke, kn) ! Output: last two entries

            end if

         end if

         ! Calculating u coefficients and source
         call get_u_coefficients(nx, ny, modvis, dt, rp, re, rn, Jp, Je, Jn &
            , ye, yk, alphae, betae, betan, gamman &
            , vle, vln, roe, ron, roa, Uce, Vcn, au)    ! Output: last entry

         call get_u_source(nx, ny, modvis, beta, dt, rp, re, rn&
            , xe, ye, xk, yk, xke, yke, xen, yen &
            , Jp, Je, Jn, roe, ron, roa, p, vle, vln &
            , Uce, Vcn, ua, u, v, cup, sup, bu)      ! Output: last three entries

         call get_east_boundary_field(nx, ny, u, ube)

         call get_east_boundary_field(nx, ny, v, vbe)

         call get_Ucbe_Vcbe(nx, ny, xe, ye, xke, yke, ube, vbe, Ucbe, Vcbe)

         call set_bcu(nx, ny, modvis, u, betan, gamman, au, bu) ! Output: last two entries

         call get_a9d_b_rescaling(nx, ny, au, bu)

         ! Calculating v coefficients and source
         call get_v_coefficients(nx, ny, coord, modvis, dt, rp, re, rn, Jp, &
            Je, Jn, xe, xk, alphae, betae, betan, gamman, vle, vln, vlp, roe, &
            ron, roa, Uce, Vcn, av) ! Output: last entry

         call get_v_source(nx, ny, coord, modvis, beta, dt, rp, re, rn &
            , xe, ye, xk, yk, xke, yke, xen, yen &
            , Jp, Je, Jn, roe, ron, roa, p, vle, vln &
            , Uce, Vcn, va, u, v, cvp, svp, bv) ! Output: last three entries

         call set_bcv(nx, ny, modvis, v, vin, betan, gamman, av, bv) ! Output: last two entries

         call get_a9d_b_rescaling(nx, ny, av, bv)

         ! Calculating SIMPLEC coefficients

         call get_internal_simplec_coefficients(nx, ny, re, rn, xe, ye, xk, &
            yk, au, av, due, dve, dun, dvn, de, dn) ! Last six are Output

         call get_boundary_simplec_coefficients(nx, ny, de, dn) ! InOutput: de, dn

         g = 1.d0/(Rg*T)

         ! Calculating the coefficients of the linear system for pressure correction

         ! g, ro, Uce and Vcn used in this subroutine are those calculated in the previous iteraction
         ! de and dn must be calculated with the coef. of the linear sytem for u and v from which
         ! u*and v*are obtained.
         call get_p_coefficients(nx, ny, dt, rp, re, rn, Jp, Uce, Vcn, ro, g, &
            de, dn, ap) ! Output: last one

         ! Solving the linear system for u

         call lu2d9(au, dl9, du9, nxy, nx, ny) ! Matrix factorization

         call fb2d9(au, dl9, du9, ny, nx, nxy, u, tolu, nitm_u, bu) ! Solving factored linear system

         call norm_l1_9d(nx, ny, u, bu, au, norm)

         ! Solving the linear system for v

         call lu2d9(av, dl9, du9, nxy, nx, ny)  ! Matrix factorization

         call fb2d9(av, dl9, du9, ny, nx, nxy, v, tolu, nitm_u, bv) ! Solving factored linear system

         call norm_l1_9d(nx, ny, v, bv, av, norm)

         ! Calculating incorrect velocities at faces (Uce*and Vcn*)

         call get_velocities_at_faces(nx, ny, dt, rp, re, rn, xe, ye, xk, yk & ! Input
         , xen, yen, xke, yke, Jp, cup, cvp, sup, svp & ! Input
         , au, av, roa, p, u, v, uea, vea, una, vna & ! Input
         , ue, ve, un, vn, Uce, Vcn)                    ! Output

         call get_Uce_Vcn_at_boundary_faces(nx, ny, ye, u, Uce, Vcn) ! Output: Uce, Vcn

         ! Calculating source of the equation for pressure correction

         ! ro, Uce and Vcn are the incorrect ones (obtained with p*);
         ! rom, Ucem and Vcnm are those of the previous iteraction
         call get_p_source(nx, ny, beta, dt, rp, re, rn, Jp, ro, ro, roa, &
            Ucem, Uce, Vcnm, Vcn, bp) ! Output: last one

         ! Mass correction cycle

         do im = 1, imax

            call set_bcp(nx, ny, ap, bp) ! Output: last two entries

            call get_a5d_b_rescaling(nx, ny, ap, bp)

            call lu2d5(ap, dl5, du5, nxy, nx, ny) ! Matrix factorization

            call fb2d5(ap, dl5, du5, ny, nx, nxy, pl, tolp, nitm_p, bp) ! Solving factored linear system

         end do

         call norm_l1_5d(nx, ny, pl, bp, ap, norm)

         call get_pressure_density_correction_with_pl(nx, ny, pl, g, ro, p) ! InOutput: last two entries

         call get_p_extrapolation_to_fictitious(nx, ny, pin, p)

         call get_density_at_nodes(nx, ny, Rg, p, T, ro)

         call get_u_v_at_real_nodes_with_pl(nx, ny, xe, ye, xk, yk, rp, pl, &
            au, av, u, v) ! InOutput: last two entries

         call get_u_extrapolation_to_fictitious(nx, ny, modvis, betan, gamman, u) ! InOutput: last two entries

         call get_v_extrapolation_to_fictitious(nx, ny, modvis, betan, gamman, v) ! InOutput: last two entries

         call get_velocities_at_internal_faces_with_pl(nx, ny, xe, ye, xk, yk, &
            due, dve, dun, dvn, pl, ue, ve, un, vn, Uce, Vcn) ! Last six are InOutput

         call get_Uce_Vcn_at_boundary_faces_with_pl(nx, ny, pl, de, Uce) ! InOutput: Last entry

         call get_density_at_faces(nx, ny, beta, ro, Uce, Vcn, roe, ron) ! Output: last two entries

         call get_T_coefficients_and_source(nx, ny, coord, modvis, beta, dt, &
            rp, re, rn, xe, ye, xk, yk, alphae, betae, betan, gamman, Jp, Je, &
            Jn, roe, ron, roa, p, pa, cp, vlp, ke, kn, Uce, Vcn, u, v, T, Ta, &
            at, bt) ! Output: last two entries

         call get_T_source_equilibrium(nx, ny, factor, app_s, Ns, rp, re, rn, &
            Jp, roe, ron, roa, Yi, hi, Uce, Vcn, bt) ! Output: last entry

         call set_bcT(nx, ny, ccTw, Tin, Twall, T, betan, gamman, at, bt) ! Output: last two entries

         call lu2d9(at, dl9, du9, nxy, nx, ny) ! Matrix factorization

         call fb2d9(at, dl9, du9, ny, nx, nxy, T, tolu, nitm_u, bt) ! Solving factored linear system

         call norm_l1_9d(nx, ny, T, bt, at, norm)

         call get_density_at_nodes(nx, ny, Rg, p, T, ro) ! ro is output

         call get_density_at_faces(nx, ny, beta, ro, Uce, Vcn, roe, ron) ! roe and ron are output

         if (it < itdt) then
            ! Calculating convergence coefficients
            maxcu = get_maxc_9d(nx, ny, au)
            maxcv = get_maxc_9d(nx, ny, av)
            maxct = get_maxc_9d(nx, ny, at)
            maxcp = get_maxc_5d(nx, ny, ap)
            maxc = max(maxcu, maxcv, maxct, maxcp)
         end if

         ! Saving backup data
         if (mod(it, wbkp) == 0 .and. optm == 0) then

            ! Finishes cpu time
            call cpu_time(tcpu2)

            call save_backup(nx, ny, it, norm1, tcpuo + tcpu2 - tcpu1, p, pin, &
               pl, T, ro, roe, ron, u, v, ue, ve, un, vn, Uce, Vcn, pmod, &
               Ns, Yi) ! Output: none

         end if

         ! Printing results to listing file
         if (it == ito) then

            if (it == 1) norm1 = norm

            ! lid - listing file id
            write(lid,*)
            write(lid,"(A7, 4(A23))") "it", "norm/norm1", "fmi (kg/s)", &
               "fme (kg/s)", "Fd (N)"
            write(lid,*)

            write(*,*)
            write(*, "(A18, A7, A13, 3(A18))") "sim_id", "it", "norm/norm1", &
               "fmi (kg/s)", "fme (kg/s)", "Fd (N)"
            write(*,*)

            ! rid = residual file id
            write(rid, "('#', A5, 4(1X, A14))") "it", "DFd*", "norm/norm1", &
               "dt", "maxc"

         end if

         ! Printing results to listing file
         if (mod(it, wlf) == 0) then

            call get_mass_flow_rate_and_thrust(nx, ny, re, roe, u, Uce, fmi, &
               fme, Fd) ! Output: last three entries

            dFd = dabs(Fd/Fd1D - Fda_a)

            write(lid,"(I7, 4(ES23.15))") it, norm/norm1, fmi, fme, Fd

            if (len(trim(adjustl(sim_id))) <= 18) then
               write(*,73) trim(adjustl(sim_id)), it, norm/norm1, fmi, fme, Fd
               73 format(A18, I7, ES13.5, 3(ES18.10))
            else
               write(*,74) trim(adjustl(sim_id)), it, norm/norm1, fmi, fme, Fd
               74 format(A, I7, ES13.5, 3(ES18.10))
            end if

            write(rid, "(I6, 4(1X, ES14.7))") it, dFd, norm/norm1, dt, maxc

            Fda_a = Fd/Fd1D

            call flush(lid)
            call flush(rid)
         end if

         if (it < itdt) call get_time_step_increment(maxc, maxcc, dt)

         ! Parada com base no valor 'não definido' da variável norm/norm1
         if (norm/norm1 .NE. norm/norm1) then

            write(*,51)
            write(lid,51)
            write(rid,51)
            flush(lid)
            flush(rid)
            close(lid)
            close(rid)

            if (optm /= 0) then

               call depp_save_fitness(0.d0, 1, "ERROR: divergence")

            end if

            stop

         endif

         if (optm /= 0 .and. it == itmax) then

            call depp_save_fitness(0.d0, 0, "ERROR: ITMAX reached")

            stop

         end if

      end do

      ! Post-processing variables: entropy and corrections at corners
      call eq_entropy_and_corners_correction(unit_au, nx, ny, cp_type, ccTw, &
         cmod, Ns, Nr, it, itmax, itimax, tol_ej, itemax, tol_n, po, T0, OF, &
         dpY, Twall, p, pl, M, cp, Rgp, hi, ej, T, gcp, gcpf, Yi, ent, entin) ! Output entries

   end subroutine get_local_equilibrium

   !============================================================================

   ! Gets the new time step in order to ensure and accelerate the convergence. Model TSI08.
   subroutine get_time_step_increment(c, cc, dt)
      implicit none
      real(8), intent(in) :: c      !< Coefficient of convergence
      real(8), intent(in) :: cc     !< Coefficient of convergence (reference value)
      real(8), intent(inout) :: dt  !< Old/New time step (s)

      ! Inner variables
      real(8) :: h
      real(8) :: mean

      mean = (cc + 1.d0)*0.5d0

      if (c < cc) then
         h = sqrt((cc - c)/cc)
      else if (c > mean .and. c < 1.1d0) then
         h = -(c - mean)**0.3d0
      else if (c >= 1.1d0) then
         h = -0.5d0
      else
         h = 0.d0
      end if

      dt = dt*(1.d0 + h)

      if (dt > 1.d-5) dt = 1.d-5

      if (dt < 1.d-7) dt = 1.d-7

   end subroutine get_time_step_increment

   !============================================================================

end program
